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Symmetry of photoexcited states with two photoinduced carriers in two-dimensional Mott 
insulators is examined by applying the numerically exact diagonalization method to finite- 
size clusters of a half-filled Hubbard model in the strong-coupling limit. The symmetry of 
minimum-energy bound state is found to be s-wave, which is different from a d x 2_ y 2 wave 
of a two-hole pair in doped Mott insulators. We demonstrate that the difference is originated 
from an exchange of fermions due to the motion of a doubly occupied site. Correspondingly 
large-shift Raman scattering across the Mott gap exhibits a minimum-energy excitation in the 
Ai (s-wave) channel. We discuss implications of the results for the Raman scattering and other 
optical experiments. 

KEYWORDS: Photoexcited state, Hubbard model, Mott insulator, Raman scattering, Exact diagonaliza- 
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1. Introduction 

The charge gap in Mott insulators is a consequence of 
strong electron correlation represented by large on-site 
Coulomb interaction. The correlation induces novel phe- 
nomena in terms of the interplay of charge and spin de- 
grees of freedom. 1 Photoexcitation across the Mott gap 
induces two carriers, an unoccupied site and a doubly 
occupied site of electrons. In two dimensions (2D), the 
motion of the two carriers is strongly affected by the pres- 
ence of localized spins in the background: The propaga- 
tion of a carrier is known to induce a spin cloud around 
the carrier as a consequence of the misaligned spins along 
the carrier- hopping paths, and the two carriers prefer to 
form bound states by minimizing the loss of magnetic 
energy due to the spin cloud. 2 The formation of such 
bound states is similar to the case of two holes introduced 
into 2D Mott insulators. The bound state formed in the 
two- hole ground state has d x 2_ y 2 symmetry (equivalently 
the Bi representation in the D4 group), provided that 
the exchange interaction between localized spins is not 
small compared with the nearest-neighbor hopping am- 
plitude. 3 This is considered to be related to d x 2_ y 2 su- 
perconducting symmetry in the high-T c cuprates. 

The nature of photoexcited states in 2D Mott insu- 
lating cuprates such as Si^CuG^C^ has been examined 
by using linear 4,5 and nonlinear optical response exper- 
iments 6-9 and large-shift Raman scattering experiments 
across the Mott gap. 10,11 Among these experimental 
techniques, the large-shift Raman scattering is the best 
one to see the symmetry of the photoexcited states di- 
rectly. However, not only excitations related to the bound 
states in the photoexcited states but also d-d transitions 
between 3d orbitals contributes to the Raman scatter- 
ing. 11 Therefore, it is important to give information on 
the contribution of the bound states to the large-shift 
Raman scattering. 

There are several theoretical studies related to the pho- 
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It has been shown analytically and numer- 
ically 14 that the bound states with odd-parity (dipole- 
allowed states) are higher in energy than those with even- 
parity (dipole- forbidden states). However, symmetry of 
the lowest-energy bound state has not been clarified yet, 
although two possibilities, the A2 12 or B1, 13 have been 
proposed. Therefore, it is important to clarify the sym- 
metry of the bound sate in order to fully understand the 
nature of the photoexcited states. Also it is interesting 
whether the symmetry is related to that of a bound state 
of two-hole pair in doped Mott insulators. 

In this paper, we theoretically examine symmetry of 
photoexcited states and large-shift Raman scattering in 
the 2D Mott insulators. We apply a numerically exact di- 
agonalization method to finite-size clusters of a half-filled 
Hubbard model in the strong-coupling limit. In the cal- 
culations, we introduce various boundary conditions. An 
averaging procedure over twisted boundary conditions 
is also used to reduce finite-size effects. The symmetry 
of the lowest-energy bound state is found to be neither 
A2 nor Bi, but Ai (s-wave symmetry). We demonstrate 
that, if sign changes due to an fermion exchange caused 
by the motion of a doubly occupied site are not taken 
into account, the symmetry becomes Bj. The large-shift 
Raman scattering exhibits a lowest-energy excitation in 
the Ai channel. This is different from the experimental 
data showing a lowest-energy excitation with A2 sym- 
metry 10, 11 and thus supports a proposal that the A2 ex- 
citation is due to a d-d transition from d x 2_ y 2 to d xy 
orbitals. 11 

The rest of this paper is organized as follows. We intro- 
duce an effective Hamiltonian of the half-filled 2D Hub- 
bard model in the strong-coupling limit, and show out- 
lines of the procedure for choosing boundary conditions 
in § 2. In § 3, calculated results of the distribution of 
photoexcited states are shown and the symmetry of the 
lowest-energy bound state is discussed. In § 4, we show 
calculated spectra of the large-shift Raman scattering 
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and compare them with experiments. The summary is 
given in § 5. 

2. Model and Method 

The Hubbard model is given by 

#Hub = H t + Hu (1) 

with 

fft = -i£(ct CT c i+5;CT +H.c.) ( 2 ) 

and 

H u = uY / ^n hl , (3) 

i 

where c\ a is the creation operator of an electron with 

spin <r at site i, rii^ — c\ a Ci t(T , the summation of 5 runs 
over a pair of the nearest-neighbor sites between i and 
i + 5, t is the nearest-neighbor hopping integral, and U is 
the on-site Coulomb interaction. In the strong coupling 
limit (U » t), there is no doubly occupied site in the half- 
filled ground state, resulting in the Heisenberg model: 

H = J^2(S 1 -S 1+S -^) , (4) 

where Si is the spin operator with S = 1/2 at site i, 
rii = n^f + n^j, and J = 4t 2 /U. 

Photoexcited states across the Mott gap have both 
one doubly occupied site and one vacant site. An effec- 
tive Hamiltonian describing the photoexcited states is 
obtained by restricting the Hilbert spaces to a subspace 
with one doubly occupied site. By performing the second 
order perturbation with respect to the hopping term H t , 
eq. (2), the effective Hamiltonian is given by 14,15 

h cS = UiHtUi - ini jr t n 2 ff t iii + ^HtUoH^ + u 

(5) 

where Ho, IIi, and H2 are projection operators onto 
the Hilbert space with zero, one, and two doubly occu- 
pied sites, respectively. We note that the effective model, 
eq. (5), can reproduce well the optical conductivity of 
the Hubbard model with large U under periodic bound- 
ary conditions for the N = 18 and 20 clusters 17 and 
under antiperiodic boundary conditions for the TV = 16 
cluster. 18 Hereafter h — e — c — 1, e and c being the ele- 
mentary charge and the speed of light, respectively, and 
the distance between the nearest-neighbor sites in the 
two-dimensional lattice is set to be unity. Throughout 
this paper, we take U/t = 10. 

We use the exact diagonalization method based on the 
Lanczos algorithm to calculate low-lying eigenstate of the 
Heisenberg and effective Hamiltonians. In 2D systems, 
one uses a iV-site square lattice with the translational 
vectors R a = lx + my and = — mx + ly, being that 
N = l 2 + m 2 with integers I, m > 0. Here, x and y are the 
vectors connecting nearest-neighbor sites in the x and y 
directions, respectively. In this study, we take N = 18 
(Z = 3, to = 3), N = 20 (I = 4,m = 2), and N = 26 
(I = 5,to= 1). 
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We impose periodic boundary conditions for the clus- 
ters along both the x and y directions. In addition, an- 
tiperiodic boundary conditions and mixed boundary con- 
ditions (one direction is periodic and the other is an- 
tiperiodic) are used to check the effect of the boundary 
condition on physical quantities. In such small-size clus- 
ters, we are not free from finite-size effects that some- 
times make the results unreliable. In order to reduce 
the finite-size effects, we also introduce various bound- 
ary conditions with twist and average physical quantities 
over the twisted boundary conditions. This procedure 
has been applied for various quantities in the t-J 19 and 
t-t'-f'-J 2 " models. 

The twist induces the condition that c i+ # ai<J = e I(?ia Cj i<T 
and c i+ R b . a = e l ^ h Ci. a with arbitrary phases <j) a and fa. 
Note that 4> a = 4>b = (it) corresponds to the periodic 
(antiperiodic) boundary conditions. The phase 4> a (b) is 
defined as 4>a(b) = K ' ^0(6) with an arbitrary momentum 
k = k x x + K y y. k usually scans an area surrounded by 
a square with four corners at (k x , n y ) — ±^(Z — to, I + 
to) and ±7f(i + to, — Z + to). In order to perform the 
averaging procedure, we choose many k in the square 
with equal intervals of 7r/45, 7r/40, and \/2tt/8^/TS for 
the N =18, 20, and 26 clusters, respectively. The total 
number of k, N k , results in N K =450, 320, and 64 for 
N =18, 20, and 26, respectively. The reason that N K 
decreases with increasing N is that time and memory for 
computing increase with N. In particular, for N = 26 we 
need to use supercomputing systems and the computing 
time of the optical conductivity for a k point is about 
130 minutes by using Hitachi SRI 1000 at the Institute for 
Solid State Physics, the University of Tokyo. Therefore, 
it is time-consuming to increase N K for N — 26. We note 
that imposing the twist is equivalent to transforming the 
operator cj ^c^s^ into exp(i« • S)c\ a c i+ s, a , 8 being the 
displacement vector from site i to i + S. Therefore, the 
Heisenberg model, eq. (4), is independent of the choice 
of K. 

3. Symmetry of Photoexcited States 

The regular part of the optical conductivity detects 
optical-allowed photoexcited states. The real part of the 
conductivity for a given k is expressed as 

where |0) is the ground state of the Heisenberg model 
with energy E , and I^J^) represents a photoexcited state 
with energy j£ is the x component of the current 
operator under the twist up to second order of t: 

2,<5,cr i,S.S' .a 

( iK -(5'-5)-t ~t ~ ~ 

-e iK ' 5 XW4-<^+<5,-<^+<5v) • ( ? ) 

The creation and annihilation operators, c\ a and c, ;Cr , 
are projected onto the subspace with either zero or one 
doubly occupied site. The 6 X is the x components of S. 



Full Paper 



J. Phys. Soc. Jpn. 



Full Paper 



N=20 



N=26 



0.3 r 



0.2 - 

b" 0.1 



. (a) 









- (b) 




J 






A/=20 



T. Tohyama 3 

N=26 



0.3 


- (a) 


■ 


0.2 






\ 1* 0.1 






0.0 







(b) 





Fig. 1. (Color online) Optical conductivity <j k (u>) for half- filled 
Hubbard clusters in the strong coupling limit with the size of 
N = 20 (a) and N = 26 (b), and dynamical correlation func- 
tion S K (w) of Ai-symmetry operator (the thick line) and of Bi- 
symmetry operator (the thin line) for N = 20 (c) and N = 26 
(d). C//t = 10. Periodic boundary conditions are employed and 
delta-functions are broadened by a Lorentzian with a width of 
0.051 



For the periodic boundary conditions, k is chosen to sat- 
isfy k ■ R a = k ■ R/, = (tt). For the mixed boundary 
conditions, we take an average of <j k (lo) with k ■ R a = 
and k ■ Rb = 7r and that with n ■ R a = tt and k ■ Rb = 0. 
The optical conductivity under the averaging procedure 
is given by 

1 V 



Fig. 2. (Color online) The same as Fig. 1, but under antiperiodic 
boundary conditions. 



N K 



(8) 



We show the optical conductivity under various 
boundary conditions for the N = 20 and N — 26 clusters 
in the panels (a) and (b), respectively, of Figs. 1-4: Re- 
sults for the periodic, antiperiodic, mixed, and averaged 
boundary conditions are shown in Fig. 1, Fig. 2, Fig. 3, 
Fig. 4, respectively. Both the N — 20 and N = 26 cases 
in Figs. 1(a) and 1(b) show a globally similar distribution 
of the spectral weight, exhibiting two prominent struc- 
tures: One is a broad-peak structure centered at around 
lo = 10<, and the other is an absorption-edge structure at 
around uj = 6t separated from the broad peak. The latter 
structure is sensitive to the exchange interaction J, 2 and 
thus originates from magnetically induced bound states. 
This global distribution of the spectral weight does not 
change under the antiperiodic and mixed boundary con- 
ditions as shown in Figs. 2 and 3, though fine structures 
in the spectra depend on the boundary conditions. 

The similarity of the spectral-weight distribution be- 
tween N = 20 and N — 26 is also seen in Figs. 4(a) and 
4(b), where we take the averaging procedure, eq. (8). 
The averaging is found to reduce a difference of the edge 
structure for N = 20 and N — 26 under the periodic 
boundary conditions 21 [see Figs. 1(a) and 1(b)]: A single 
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Fig. 3. (Color online) The same as Fig. 1, but under mixed 
boundary conditions, where the average of two cases [periodic 
(antiperiodic) along the x (y) direction and antiperiodic (peri- 
odic) along the x (y) direction] is taken. 



peak appears at to = 6.2t for TV = 20, while for N = 26 
there are four peaks with similar magnitude within the 
region of 6t < lo < 7t. After taking the averaging over n, 
the edge structure becomes wider for N = 20 but nar- 
rower for N = 26. Therefore, the difference appearing in 
the cases of the periodic boundary conditions is reduced 
in Fig. 4, although a two-peak structure still remains for 
N = 26. The reduction demonstrates efficiency of the av- 
eraging procedure. We note that cr avc (uj) of TV = 18 (not 
shown) is similar to that of N — 20. 

In order to examine whether optical-forbidden states 
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Fig. 4. (Color online) The same as Fig. 1, but an averaging pro- 
cedure for twisted boundary conditions is employed. 



are lower or higher in energy than allowed ones, we in- 
troduce operators with Ai and Bi symmetry composed 
of the nearest-neighbor hoppings: 

(9) 

where the plus (+) and minus (— ) signs correspond to 
Ai and B 1; respectively, and IS the a component of 
k. The dynamical correlation function of the operators 
for a given k is given by 

S "H = JfY, K*m Q )\ 2 ^ E ™ + E <>) > ( 10 ) 

m 

and an averaged correlation is expressed as 

SavcM = ^]Ts K H. (ii) 

We note that the final states |^J^) cannot be divided into 
subspaces of these representations because of the twisted 
boundary conditions. 

S K (w) for N = 20 and N — 26 is shown in the pan- 
els (c) and (d), respectively, of Figs. 1-3, and S ave (u}) 
is shown in Figs. 4(c) and 4(d). For the cases of the 
periodic, antiperiodic, and mixed boundary conditions 
(Figs. 1 - 3), we find edge-structures in the Ai correla- 
tion separated from broad spectrum, indicating the pres- 
ence of bound states. We also find that in these boundary 
conditions the energy of the edge position of the Ai cor- 
relation is lower than those of the Bi correlation and 
u K (ui). As is the case of o" ave (w), the averaging procedure 
makes differences between N = 20 and N = 26 smaller. 
After averaging, the fact that the edge of Ai is lower than 
that of Bi is still preserved, indicating an intrinsic na- 
ture of the photoexcited states of the half-filled Hubbard 
model with large U. As will be shown below, the states 
having A2 and B2 components are higher in energy than 




Fig. 5. (Color online) Dynamical correlation function Sav C (w) of 
Ai-symmetry operator (the thick line) and of Bi-symmetry op- 
erator (the thin line) for a TV = 20 half-filled Hubbard cluster 
in the strong coupling limit, (a) U/t = 10 and V/t = 1. (b) 
U/t = 10, but sign changes due to the exchange of fermions 
caused by the motion of a doubly occupied site are neglected. 
An averaging procedure for twisted boundary conditions is em- 
ployed and delta-functions arc broadened by a Lorentzian with 
a width of 0.05t. 



the Ai dominated states. These results indicate that the 
lowest-energy bound state in the photoexcited states of 
the Hubbard model has the Ai symmetry, i.e., s-wave 
symmetry. 

This conclusion is not altered even if we include an 
attractive interaction, V, between neighboring doubly 
occupied and vacant sites, which induces an excitonic 
bound state. The resulting effective Hamiltonian is given 

by 

H v = H eS - V \] [Vit^Tfyj, (1 - 7tj + a,-|-) (1 - n i+ s,i) 

1,8 

+ (1 - n i.{) (1 - n i,i) n i+6,l n i+S,l] ■ ( 12 ) 

Figure 5(a) exhibits S ave (w) for the N = 20 cluster with 
V = t. We find no qualitative change of the bound states 
as compared with Fig. 4(c): The excitation energy of the 
Ai bound state is lower than that of Bi. 

In the two-hole doped Mott insulator, the hole pair 
forms a d x i_ y i wave in the ground state. 22 This is in 
contrast with the present results that the two-carrier 
pair produced by photoexcitation forms an s wave. It 
is important to clarify what is the origin of this dif- 
ference. We can easily notice a remarkably difference 
between the two-hole and two-photoexcited-carrier sys- 
tems. That is the difference of the electric charge of 
carriers. In the photoexcited states, one of the carriers 
is not a hole but contains two electrons. Therefore, in 
a basis representation where electrons are sorted in or- 
der of site, the motion of this carrier inevitably induces 
an exchange of fermions and gives an extra sign. For 
instance, in the case of a two-site cluster with a dou- 
bly occupied site and a singly occupied site, the hop- 
ping of the doubly occupied site induces an extra sign: 
(— tc\ ^ci^)c\ ,c\ ,cj 1 I vac) = +tc| ic| »c| , |vac), |vac) 
being the vacuum state. 

In order to check the effect of the fermion exchange, 
we introduce an effective Hamiltonian that is the same 
as eq. (5) but reversing the sign of hopping of the doubly 
occupied site. Figure 5(b) shows 5 ave (w) for the effective 
model. We find that the Bi correlation exhibits bound 
states whose energies are lower than those of Ai, indi- 
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Fig. 6. (Color online) Large-shift Raman scattering R K (u>) for a 
N = 20 half- filled Hubbard cluster in the strong coupling limit. 
U/t = 10. (a) Ai, (b) A 2 , (c) Bi, and (d) B 2 scattering. The 
thick and thin lines represent the results under periodic and an- 
tiperiodic boundary conditions, respectively. Delta- functions are 
broadened by a Lorentzian with a width of 0.05t. 



eating the lowest-energy state to be d x 2_ y 2-w&ve. From 
this result, we can conclude that the fermion-exchangc 
process for the doubly occupied site plays a crucial role 
in making the Ai bound state lowest in energy. 

4. Large-Shift Raman Scattering 

Large-shift Raman scattering for the half-filled Hub- 
bard model is a good quantity to study the photoexcited 
states in the symmetry-resolved form. The Raman inten- 
sity for a given k is given by 



f 



(13) 



with symmetry-resolved Raman operators Mr: Mai 



M xx - 

Mb2 — ^xy 

expressed as 



M vv , M A2 = M 



L yyi 

M. 



M x 



xy — mp, Msi 

The operator M a p (a, (3 



M, 



a/3 



J f3 l*r 



\3S 



En 



"7 



M M , and 

= a;, y) is 



(14) 



where Jg connects two subspaces in the photoexcited 

states and is given by J" = j% + i[U~ 1 HiH t n.oH t n.i, 0], 

where is the component of the total position opera- 
tor. The Raman intensity under the averaging procedure 
is given by 



(15) 



We set the incident photon energy u>i to be u>i = U+6t = 
16t and focus on the Raman shift u> of order of U. We 
take the relaxation energy rj to be 0.2t. 

Figure 6 shows R K (u>) for the N — 20 cluster with the 
periodic and antiperiodic boundary conditions. Although 
global features up to the high-energy region are similar 
between the two boundary conditions, edge structures 
of the spectrum depend on the boundary conditions. In 



Fig. 7. The same as Fig. 6 but an averaging procedure for twisted 
boundary conditions is employed. 



order to reduce the dependence, we use the averaging 
procedure, eq. (15). R ave (w) is shown in Fig. 7. Among 
all possible symmetry, the Ai Raman scattering exhibits 
the lowest-energy excitation, being consistent with the 
Ai correlation function discussed above. For the Bi and 
B2 Raman scattering, spectral weights appear slightly 
above the Ai scattering. In the A2 scattering, no weight 
is observed at the edge region around w = 6t. These 
results demonstrate again that the Ai state is the lowest- 
energy state in the photoexcited states of the 2D Mott 
insulator. 

Since the N = 20 cluster is not a simple square lattice 
but a tilted one, it does not have the D4 point group. 
In order to check whether the data in Figs. 6 and 7 are 
artifacts of tilted lattice or not, we examine R K (u>) and 
-Rave(w) for the N = 18 cluster that has the D4 group. 
The results are shown in Figs. 8 and 9. In the case of the 
periodic boundary conditions where the D4 point group 
is fully held, R K (ui) exhibits the same behaviors as the 
case of N = 20 under the periodic boundary conditions: 
The Ai scattering shows the lowest-energy excitation and 
the A 2 scattering has no weight at around u> = 6t. Al- 
though the difference of R k (uj) under the periodic and 
antiperiodic boundary conditions is larger in the case of 
N = 18 than that in the case of N = 20 shown in Fig. 6, 
i?ave(w) obtained after the averaging in Fig. 9 shows spec- 
tral features similar to the case of N = 20 in Fig. 7. This 
indicates that the behaviors of Raman intensities shown 
in Figs. 6 and 7 are inherent in the square-lattice Hub- 
bard model. 

The large-shift Raman scattering experiments for 
parent compounds of high-T c superconductors such as 
Nd2CuC>4 have shown that the A2 scattering has the 
largest spectral weight whose excitation energy is lower 
than the absorption-peak position. 11 This is completely 
different from the present results shown in Figs. 7 and 
9. We consider that, as discussed in ref. 11, the observed 
large A2 scattering would come from a d-d transition 
from d x 2_ y 2 to d xy orbitals through photoexcited states. 
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Fig. 9. The same as Fig. 6 but for TV = 18 and an averaging 
procedure for twisted boundary conditions is employed. 



5. Summary 

Before summarizing our results, wc discuss implica- 
tions of the present results for recent two-photon absorp- 
tion (TPA) experiment on NCI2CUO4. 9 In the experiment, 
there is a small absorption hump whose energy is about 
0.3 eV lower than the optical (linear) absorption energy. 
Combining this experimental data with third-harmonic 
generation data 8 ' 23 and analyzing the data by using a set 
of discrete excitation-energy levels, Kishida et al. have ar- 
gued that the hump structure in TPA comes from a pho- 
toexcited state with Bx symmetry and thus the lowest- 
energy photoexcitation is due to the 13i symmetry. 23 This 
is different from the present results that the minimum ex- 
citation comes from the Ai symmetry. According to the 
Kishida's analysis, it is difficult to assign the hump as 
the Ai excitation. Therefore, a possible interpretation is 
that the hump structure in TPA would originate from a 
d-d excitation between, for instance, d x 2_ y 2 and d 3z 2_ r 2 
orbitals. However, in order to come to final conclusions, 



we need further theoretical studies of, for instance, TPA 
spectra under the presence of various d orbitals. 

In summary, we have examined symmetry of photoex- 
cited states and symmetry-resolved large-shift Raman 
scattering in the 2D Mott insulators by using an effec- 
tive Hamiltonian of a half-filled Hubbard model in the 
strong-coupling limit and a numerically exact diagonal- 
ization method on finite-size clusters. The symmetry of 
the lowest-energy bound state is found to be neither 
A2 nor Bi, but Ai. Therefore, the symmetry is differ- 
ent from that of a hole pair in doped Mott insulators. 
We demonstrate that the difference is originated from 
fermion exchange induced by the motion of a doubly oc- 
cupied site. The large-shift Raman scattering exhibits a 
lowest-energy excitation in the Ai channel. This is differ- 
ent from the experiments showing a lowest-energy exci- 
tation with A2 symmetry, which is probably due to a d-d 
transition from d x 2_ y 2 to d xy orbitals. It seems to be nec- 
essary to include other d orbitals into the photoexcited 
states for a complete understanding of the experimental 
data. This remains as a future problem. 
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